clearvars -except do_QC

date     = SAT_SST_func_IO('date');
if do_QC == 0
    P.app      = 'qcu';
else
    P.app      = 'qcf';
end

% *************************************************************************
% Read masks
% *************************************************************************
dir      = SAT_SST_func_IO('GHCN');
mask     = ncread([dir,'land_percent2_qd.nc'],'lsm');
mask     = mask([721:1440 1:720],end:-1:1);

is_coastal = zeros(1440,720);
is_island  = zeros(1440,720);

for ct_lon = 1:1440
    for ct_lat = 2:719
        
        l_lon = [-1 0 1] + ct_lon;
        l_lon(l_lon > 1440) = l_lon(l_lon > 1440) - 1440;
        l_lon(l_lon < 1) = l_lon(l_lon < 1) + 1440;
        
        l_lat = [-1 0 1] + ct_lat;
        
        temp = mask(l_lon,l_lat);
        
        if nanmean(temp([1:4 6:9])) < 10
            is_island(ct_lon,ct_lat) = 1;
        end
        
        if min(temp([2 4 6 8])) < 1 || temp(5) < 50
            is_coastal(ct_lon,ct_lat) = 1;
        end
    end
end

% *************************************************************************
% Read station metadata
% *************************************************************************
dir      = SAT_SST_func_IO('GHCN');
file_inv = [dir,'ghcnm.tavg.v4.0.1.',date,'.',P.app,'.inv'];
fid = fopen(file_inv,'r');
ch = fscanf(fid,'%c',10000000000000);
fclose(fid);

clear('inv')
list = [0 find(ch == 10)];
for ct = 1:numel(list)-1
    temp = ch(list(ct)+1:list(ct+1));
    inv(ct,:) = temp(1:68);
end
inv_ID = inv(:,1:11);
lat = str2num(inv(:,13:20));
lon = str2num(inv(:,22:30));  lon(lon<0) = lon(lon<0) + 360;
lev = str2num(inv(:,32:37));
lev(lev > 9900) = nan;

% *************************************************************************
% Read distance to the coast
% *************************************************************************
dir       = SAT_SST_func_IO('GHCN');
Table     = readtable([dir,'GHCN_Metadata_MicroSite.csv']);
names     = cell2mat(table2cell(Table(:,1)));
distance  = cell2mat(table2cell(Table(:,7)));
lon_meta  = cell2mat(table2cell(Table(:,3)));
lat_meta  = cell2mat(table2cell(Table(:,4)));
[~,place] = ismember(inv_ID,names,'rows');
dis = nan(size(place,1),1);
for ct = 1:size(place,1)
    if place(ct) > 0
        temp = distance(place(ct));
        if ~isempty(temp)
            dis(ct) = temp;
        end
    end
end

% *************************************************************************
% Exclude stations following Cowtan et al. (2018)
% *************************************************************************
l_exclude = abs(lat) > 60 | (lon > 0 & lon < 52 & lat > 28 & lat < 90) | dis <= -10 | isnan(dis);

use_station     = ~l_exclude;
island_station  = grd2pnt(lon,lat,is_island,0.25,0.25)  & ~l_exclude;
coastal_station = grd2pnt(lon,lat,is_coastal,0.25,0.25) & ~l_exclude;

ID_use_station     = inv_ID(use_station,:);
ID_island_station  = inv_ID(island_station,:);
ID_coastal_station = inv_ID(coastal_station,:);
dis = dis(use_station);

% *************************************************************************
% Save station list
% *************************************************************************
dir_save  = SAT_SST_func_IO('GHCN');
file_save = [dir_save,'Coastal_station_list_',P.app,'_',date,'.mat'];
save(file_save,'ID_use_station','ID_island_station','ID_coastal_station','dis','-v7.3')
